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ABSTRACT 


This paper presents the application of a generalized optimality criteria to framed structures. The 
optimality conditions, Lagrangian multipliers, resizing algorithm, and scaling procedures are all represented 
as a function of the objective and constraint functions along with their respective gradients. The optimization 
of two plane frames under multiple loading conditions subject to stress, displacement, generalized stiffness, 
and side constraints is presented. These results are compared to those found by optimizing the frames using 
a nonlinear mathematical programming technique. 

INTRODUCTION 

Weight optimization of large aerospace structures requires the use of efficient optimization methods due 
to the potentially excessive number of design variables and related constraints. In order for optimization to 
be seriously used by the preliminary designer, the method must be able to handle multidisciplinary problems 
(thousands of design variables and their related constraints) and must be efficient (produce designs in hours 
not days). 

In the late sixties and early seventies the optimality criteria approach to structural optimization was 
developed [l]. At that time, and in subsequent work, the optimality criterion was derived for an individual 
problem and although very efficient, it was criticized for its lack of generality. Recently, Venkayya[2) has 
generalized the optimality criteria to apply to any structural optimization problem to which the sensitivities 
of the objective and constraint functions can be computed. First, this paper will briefly state the optimality 
conditions. Next, detailed descriptions of the Lagrangian multipliers, scaling formulation along with a 
redesign procedure for stress, displacement, and generalized stiffness constraints applied to plane frames 
will be presented. Finally, the optimization results of some plane frames using the generalized optimality 
criterion will be given along with a comparison to results found by nonlinear mathematical programming. 

PROBLEM STATEMENT 

The optimization of a structure for minimum weight can be stated mathematically as: 

Minimize the objective function 


W(A) = '*Tp i l i A i (1) 

t-1 


subject to the constraints 


Z j (A)<Z j j = 1,2, ...n (2) 

I 

A < * < A^ ( 3 ) 

where W (A) is the weight of the structure, pi , A*, and are the specific weight, the cross-sectional area, 
and the length of the ith element respectively. The Zj(A) consist of all n of the behavioral constraints and 
Z 0 is the given allowable for Z 3 (A). The summation in equation (1) is over all m elements in the structure. 
However, this does not imply that all the elements are required to participate in the design algorithm. In 
addition to the constraints Z 3 , each design variable (A*) has upper bounds (A^) and lower bounds (A^) 
referred to as side constraints. 


DESIGN VARIABLES 

The design variables chosen for the plane frame are A*, U z , and S**, where A* is the cross-sectional 
area, /** is the moment of inertia about the 2 -axis and S iz is the section modulus about the 2 -axis. These 
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variables are not independent; therefore, A x is chosen as the primary variable with I xz and Si z expressed as 
explicit nonlinear functions of Ai in the form 


Iiz = aiA”* (4) 

Siz^nA? (5) 

Depending on the type of cross-section and the assumptions being made n x will vary from 1 to 3 , and v x will 
vary from 1 to 2. It is important to note that this method is general and any design variable can be chosen 
such as width of the section, thickness of the flange, or thickness of the web of the section. This section will 
focus on relations between A { , Ii z , and Si z for solid rectangular cross-sections and a three spar box section. 
Three separate cases will be presented for each type of section. Each case will make varying assumptions 
about the width of the section, depth of the section, thickness of the flange, thickness of the web, and ratios 
of these quantities. 

For the rectangular section in figure 1 having depth d and width b the following cases for the relations 
between Ai, U z , and S{ z are presented. 

Case 1: Assume b to be constant and d is allowed to vary. 


I\z — n i — 3 &i ~~ 22^2 

S\z = 7» * Vi — 2 7* = ^ 


Case 2: Assume the ratio b/d is equal to some constant C 0 and let b and d vary. 


( 6 ) 


I iz = a. A? rii = 2 
S i , = 'nA v i ‘ Vi = 3/2 


a. = 

7t = 


1 


12C a 

1 


6-y/Co 


Case 3: d is assumed to be constant while b is allowed to vary. 


Iiz — 1 aj 

Siz = l.A? f, = 1 ~ii = 


<P_ 

12 

d 

6 


( 7 ) 


( 8 ) 


Next, consider the three spar box section in figure 2. The area A*, moment of inertia I xz , and the section 
modulus can be expressed as: 


or 


or 


*=“{ 2 <i + T>^ + s ('- 7 >t} 


A x — bdCi 

\ / ^ £/ \ 3 




/ - bd3 c 

h * - l2° 2 


s '- = x{(r7b )l(1 + ‘f ,(1 + 7) s+ 


t, Id 


( 9 ) 


(10) 
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or 


Si. - ^C 2 C S 


( 11 ) 


Following these definitions the three cases can now be presented. Casel: The t w , b, and tj/d are all held 
constant. 



Ai = bdC\ 

Iiz = 

rii = 3 oti = 

= liA”‘ 

^ = 2 7t = 


C 2 


12fc 2 C? 

C 2 

6 CfC s b 


Case 2: b/d a constant C\, t w /d, and t//d a constant. 

Ai = bdCi 

Iiz = aiA^ rii = 2 at = 


i2C?c t 
C 2 c< 


Siz-nA*' v* — 3/2 Tf i Q C ^ Cl c t )s / 2 

Case 3: b allowed to vary with d, tu,/b, and t//d held constant. 

A, = bdCi 


Iiz — ’ Mi — 2 OCi — 

Siz = 1>AV Vi = 3/2 7 . = 


d 2 C 2 

12Cj 

dC 2 


( 12 ) 


(13) 


(14) 


e>c x c 3 

These relations for the three spar box can easily be extended to a n spar box cross-section. Also, the I-section 
can be show to be a special case of the spar box with n = 1. 

CONDITIONS OF OPTIMALITY 


The optimisation problem defined by equation (1) can be restated in Lagrangian form as 


LU,\) = W(A)-'£x i lZ J -Z J ) 

3 = 1 


p = number of active constraints (15) 

L is the Lagrangian function, and are the Lagrangian multipliers corresponding to the active constraints. 
A constraint will be defined as active if Zj = Z, Minimization of the Lagrangian function with respect to 
the design variable Ai gives 


dL = dW^_Y' x d _h =0 

dAi dAi J dAi 


i - 1, 2 ...m 


(16) 


Equation (16) can be rewritten as 


p 

^ ^ Cgj Xj — 1 i 1, 2..., m 

3 = 1 


(17) 
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where is the ratio of constraint to objective function sensitivities and is given by 


_ dZj/dAi 
e ° " dW/dAi 

or in matrix form 

eX = 1 

where e is an m x p matrix, A is an p x 1 matrix and 1 is a m x 1 matrix. 

LAGRANGIAN MULTIPLIERS 

The Lagrangian formulation introduces more unknowns in addition to the m design variables. These 
additional unknowns are the Lagrangian multipliers and there are as many Lagrangian multipliers as active 
constraints p. Thus, it is necessary to solve for m + p unknowns. This section discusses some methods of 
solution for the Lagrangian multipliers. 

Premultiplying equation (19) by t l w yields 

t l we X — e l ty\ ~ Z (20) 

where the weighting matrix w is an m x m positive definite diagonal matrix. Here the diagonal elements of 
the w matrix are taken to be the individual weights of each element (u>^ = piliAi) . Equation (20) can now 
be stated as 


( 19 ) 


HX = Z (21) 

Although the H matrix is non-singular, it is an implicit function of the final design variables. Thus, equation 
(21) represents a non-linear set of equations. Thus, some approximate methods are used instead of trying 
to actually solve for the As by some iterative scheme. The first approximate method considered consists of 
using the information at the current design and solving the equation (21) for the Lagrangian multipliers by 
inverting the H matrix. 

A - (22) 

One of the drawbacks of this method is that there is no guarantee that all the Lagrangian multipliers will 
have the appropriate sign. This causes problems when attempting to resize the design variables. The second 
method considered was originally developed in 1973 by Venkayya and coworkers[3j. This method assumes 
that only one constraint is active at any one time. For a multi-constraint problem the A’s simply become 
weighting parameters. For a single constraint equation (21) reduces to 


In the case of multiple constraints 


where Zj is found to be 


A = = 

(23) 

Z 


W 

(24) 

A i = =- 

Zi 


m 


= Y, ei > w " 

(25) 


»=1 


Calculating the A’s with this approximate method is computationally very efficient compared to inverting 
the H matrix. This is because wa are the weights of the individual elements and Zj is a simple function of 
the given constraint. When finding the As by this approximate method it is always assured that the As will 
have the necessary sign. 
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RESIZING ALGORITHM 


Using the optimality criterion described in equation (17) an iterative resizing algorithm can be found 
by multiplying equation (17) by A\ and solving for Ai 


a; +i 


a: 


y! 


j = i 




( 26 ) 


where £ is defined as a step size parameter and r indicates the rth iteration. For most problems £ = 2 is 
chosen and gives a good rate of convergence. 

SCALING PROCEDURE 


Once a new design has been generated by the resizing algorithm the constraint surface must be found. 
This is done by the use of a scaling procedure. Let A be the current design vector with the new design found 
by A = Ai4 where A is the scale factor. If dA is the difference between two designs, it can be written as 

dA = A~A = (A -1)4 (27) 

If the response of the structure is R , then performing a first order Taylor Series expansion on R about the 
current design point A yields 

R = R + f^—dAi (28) 

t=l ° ' 

Substituting equation (27) for dA y dividing both sides of the equation by the response R , and realizing that 
R ~ R = dR yields 


dR 

IT 


(A-l) 


E m 

i = 1 


SR A. 
SA t 


R 


(29) 


In equation (29) the term / R is either < 0 or > 0 depending on the type of constraint being 

considered. For this work only stress and displacement constraints are being investigated and for these cases 
[jZi^i /& * s — 0. Now defining // as 



\^m dR a 
= l dA % 

R 

(30) 

equation (29) can be expressed as 

^ = (i-a)m 

(31) 

Solving equation (31) for the inverse 

of the scale factor A and defining b — ^ —■ gives 



I - 1 6<k 1 

A 1-6 

(32) 

By performing a binomial expansion 

and ignoring higher order terms equation (32) becomes 



J — ‘ 

(33) 

Rearranging equation (33), adding 1 to both sides, and defining as the target response ratio equation (33) 

becomes 


(34) 
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where fi is 


New Response R + dR 
Initial Response R 

Finally, solving equation (34) for the scale factor gives 


( 35 ) 


A=-^— 

P + M “ 1 

In the case of truss or membrane structures fi = 1 and A reduces to lj fi which is the exact scale factor for 
stress and displacement constraints. 

SPECIALIZATION TO BENDING ELEMENTS 


In the following sections the A j, and A for bending elements subject to stress, displacement, and 
generalized stiffness constraints will be discussed. 

DISPLACEMENT CONSTRAINTS 


To find the e X j for displacement constraints the gradient of the constraint with respect to the design 
variable (dZj/dAi) is required. There are several methods for finding these gradients (finite difference, 
direct differentiation, virtual load method), but for this work the virtual load method was incorporated. 
The virtual load method consists of expressing the active constraint Zj in terms of a virtual load vector F* 
and the global displacements u. Thus, the displacement constraint Zj = Uj can be written as 

Zj = Uj = F\u (37) 


where Fj is the virtual load vector in which Fj — 1 for i = j and Fj = 0 if i / j. The e x j can be found 
by first partitioning the element stiffness matrix Ki into axial K A % and bending Kb% components. Next, 
substitute the relation a x A for 7 iz and note that dW /dA x — pi U. Then e x j for a displacement constraint 
becomes 


Pi It 


(38) 


where f\ is the virtual displacement vector corresponding to the virtual load vector F\ and is obtained from 

- 3 ~ J 

the relation 




“L 


(39) 


Once e x j is known the Lagrangian multiplier A j for the j th active constraint can be found by using 
equation (23). The resulting A’s are 

W 

A, = - =— ; — r (40) 


where the parameters pt A j and fiBj are 


Zj(p A j -H ftBj ) 

E, m =1 


tAj = 


HBj 




F)u 


(41) 

(42) 


For the scaling factor of equation (36) the parameter fi can be broken up into axial and bending parts 
HAj, HBj- Therefore, 

M = HAi + (43) 
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and the scale factor for displacement constraints becomes 


Pa j + PBj 
ft + PAj + PBj ~ 1 


(44) 


1//3 which is the scale factor for membrane structures originally found by Venkayyain 197l[l]. By inspection 
of fiAj and nsj and recalling the limits on rti for bending elements, it follows that 

1 < PAj + PBj < 3 (45) 


STRESS CONSTRAINTS 


For bending elements the stress in the j th member is expressed in terms of its bending and axial 
components as; Oj — 4 + <t 3 b- The gradients of the stress constraints (dZj/dAi) are found by the adjoint 

variable method which is a generalization of the virtual load method. That is, the constraints are recast in 
terms of a virtual load vector Fj and the global displacement vector u. The stress in a given member was 
written by Venkayya [2] as 

= TiQi (46) 

where the vector T ■ is defined as 


TJ = 


SGN n SGN n n n 
—7— 0 — — 0 0 0 
. Aj Sj 


end A 


T‘ = 


SGN „ SGN 
0 0 0 — - — 0 




Sj 


end B 


(47) 

(48) 


where SGN is the sign on the entries of the element force vector Q. and Sj is the section modulus defined 

as Sj = ")j A V j j . The element force matrix can be expressed in terms of the global displacements u, the local 
element stiffness k- , and a transformation matrix a- as 


Qi = fcj “ j « ( 49 ) 

Now the stress Oj can be written in terms of a virtual load vector Fj and the global displacements as 

<Tj = F)u (50) 

where the virtual load vector F * — T^k-a^ Unlike displacement constraints the derivative of the virtual load 
vector with respect to the design variable is not equal to zero (dFj/dAi / 0). This fact causes somewhat 
more complicated expressions for A j and A. For stress constraints becomes 


tiiWjQi -fjQj-QKAi + niKste 


where 6{j is the Kronecker delta and T 3 and Q are defined as 


ll 

1 

Q>| 

(52) 

Qj = + n ih B ))ij u i 

(53) 

Now, the Lagrangian multiplier for the ;th constraint can be found as 


W 

(54) 

\j — - — 

Zj(PAj + PBj ~ Pj) 
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where tiAjiPBj are the same as those for the displacement constraints. pj is a new term introduced due to 
the fact that dF^/dAi ^ 0 for stress constraints and it is found to be 


H = 



(55) 


Finally, following the derivation of equation (36) the scale factor for bending elements subject to stress 
constraints can be found to be 


JMj H- fiBj ~ 

P + PAj + PBj ” “ 1 


(56) 


GENERALIZED STIFFNESS CONSTRAINT 


If P is the generalized forces and u is the corresponding generalized displacements then the generalized 
stiffness constraint can be written as 

Zj(A) = -Pftii i = 1, 2, ...load cases (57) 

2 

The Cij y HAj> f*Bj> and A j can be found to be 


iu'jKAi + niKBilu 

Pi U Ai 


_ utKAili 

^ - h 

E tliv} KqxU 

"PUT 

t=l 3 


Zj(tlAj + PBj) 

It is worthwhile to note that the generalized stiffness constraint does not require the use of the virtual load 
and displacement vector. This is because the information needed for the gradient of the constraint is already 
available and no new computations are needed. 

MEETING THE CONDITIONS 
OF OPTIMALITY 


The conditions of optimality state that the product between e tJ \j summed over all active constraints 
should be equal to unity at the optimum design. The are the ratios of the constraint gradients to objective 
gradients. These gradients are taken with respect to each active design variable A { . A design variable is 
considered active if it satisfies the following criteria: 


1. The variable is chosen to participate in the design iteration. 

2. The variable is within the given allowable limits. 

3. The sensitivity of the jth active constrain with respect to the design variable Ai is negative (|^~ < 0). 

It is important to note that this third criteria is constraint dependent. 

If the design variable does not satisfy the above criteria it is considered passive and the conditions of 
optimality will not be satisfied for that particular design variable. That is ma y no ^ ec l ua l to 

unity at the optimum design. If the variable does not satisfy the first two criteria then the variable is simply 
eliminated from the design for that particular iteration. However, if the variable passes the first two criteria 
there is some question on how to handle the third criteria since it is constraint dependent. It is easy to 
see that Ai may be passive for one constraint ( wtl > 0) but active for another constraint ( 1 < 0). 
This raises the question of how does one deal with the tij when a member is passive relative to a particular 
constraint? In this work if dZjdf Ai > 0 then for that particular constraint and design variable was set 
to zero. 
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MODIFYING THE RESIZING ALGORITHM 


When using the resizing algorithm equation (26) as it appears (that is taking the sum of all the CtjAj) 
it was found that this tends to over constrain the problem. The converged optimum was well above the 
known optimum. This was particularly true for multiple loading conditions. Here, instead of using the 
entire sum for resizing the maximum value of tijXj for the particular design variable was chosen. This can 
be interpreted as each variable being resized based on the constraint that is most critical for that particular 
element. This method was found to work well and allowed the algorithm to converge to the known optimum. 

REDUCING THE ERROR 
IN THE SCALE FACTORS 

Due to the Taylor Series and binomial approximations the scaling factors in equations (40,56) are only 
valid within certain limits. This is especially true when the structure is primarily in bending. This is because 
equations (40,56) do not reduce to the exact bending scale factor (l//?) 1 ^ if axial contributions are ignored. 
It is desired that the limits on extend indefinitely without allowing the error in the scale factor or the 
response to exceed 5%. If this can be accomplished, then no additional detailed analyses are required to 
scale the design to the constraint surface. Venkayya[3] achieved this by writing an interaction formula in the 
non-dimensional parameter space /i*. Since the limits on fisj and A are known this is easily accomplish. 
The and indicate what portion of each scale factor (A axiah X bending ) must be used to generate the 
scale factor for the combined axial bending case. A linear interpolation was used and the error on the scale 
factor and response was found to be < 2% regardless of the /?. The scale factor can be represented by a 
linear interaction formula as 



where fiAj are the non-dimensional parameters found in equations (41,42), and Ji A ■ = 1 and = n. 
One could also fit a higher order polynomial between the two ranges of A to totally eliminate the error, but 
an error < 2% is generally sufficient. If each element in the structure has a different n then equation (59) 
can be written as ^ 

1 _ ( 1 \ , MB; ( 1 \ g (an , 

where now Ji B = n and n = with 

3 


SAMPLE PROBLEMS AND RESULTS 

On the basis of the preceeding derivations, a computer program written in FORTRAN 77 was developed. 
The frames used in these examples are steel structures with a specific weight of .283 lb/cubic inch and a 
modulus of elasticity of 29ksi. The type of section used is I-sections and the values for <*j, n*, ^ and 
are .2072, 3.0, .393, 2.0 respectively for all members. All problems were solved on a VAX 8600. For these 
problems the resizing was based on the generalized stiffness and displacement constraints where the scaling 
was done with respect to the stress and displacement constraints. 


^Aj \PJ V B j V / 9 


£, m 


VBj 


n =i /;*. 

F)u 


Bil 




(61) 


The first example optimized is a ten-story symmetric frame show in figure 3 that was reported by Tabak 
and Wright [4] . In this work the distributed loads used in reference [4 J were replaced by concentrated load 
(figure 4), thus creating new nodes at the mid span of each floor. By doing this the thirty member plane 
frame reported in reference (4) became a forty member structure. This frame was optimized with stress 
constraints of 22ksi on each element and displacement constraints of two inches in the horizontal direction 
for all the nodes. Figure 5 shows that the math programming method converged [5], [6], [7] in seven iterations 
to a final weight of 35,051 pounds in 37.28 cpu seconds while the optimality criteria converged to 36,421 
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pounds in four iterations with a cpu time of 4.33 seconds [8]. The optimality criteria final weight is slightly 
higher (4%) but the cpu time is significantly less (over eight times less) than that of math programming. 

The final example optimized is the 313 member frame in figure 6. This frame is subject to five loading 
conditions (figure 7), along with stress and displacement constraints. The displacement constraints are 4.0 
inches in the vertical and 12.0 in the horizontal direction at all nodes. The limit for the stresses in each 
element was 29 ksi. In figure 8 it can be seen that math programming converges to a final weight of 120,419 
pounds in fourteen iterations with a cpu time of 58 minutes. Where the optimality criteria converged to a 
final weight of 125,166 pounds in twenty-five iterations using approximately 8 minutes of cpu time. Again, 
the optimality criteria converges to a slightly higher weight (4% higher). The optimality criteria took a 
significantly large number of iterations to converge compared to math programming. This poor convergence 
is partly due to constraint switching. Even with the large number of iterations the cpu time for the optimality 
criteria algorithm is much lower than that for math programming. 

CONCLUSIONS 

The generalized optimality criteria presented in this paper can be applied to any structural optimization 
problem and related constraints provided that the constraints and their respective gradients are available. 
The math programming method finds a new design by adding and subtracting gradient information to the 
current design. Searching from point to point can be a very long and costly procedure. On the other 
hand, in the optimality criteria a redesign is computed by multiplying (not adding) gradient information 
to the current design thus, sweeping the design space instead of performing a point search. The optimality 
criteria is also fairly independent of the the number of design variables, thus allowing literally thousands 
of independent design variables. This is not the case with math programming where an upper limit on 
the number of independent design variables is between three to four hundred. When this limit is exceeded 
computer time becomes excessive and convergence is uncertain. For these problems the math programming 
method although computationally heavy gave a smooth rate of convergence and overall very good results. 

There are some disadvantages to using the optimality criteria which are evident in the sample problems. 
The optimality criteria, with the current implementation, converges to a 2% to 4% higher weight than that 
found by the mathematical programming method. Also, for the 313 member frame the optimality criteria 
method took a very large number of iterations to converge to the optimum. This is very undesirable since 
detailed analyses for large problems can become extremely expensive. 
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Load case 1 


Node 

x-Load(lbs) 

y-Load(lbs) 

1,3 


-6000. 

2, 


-24000. 

4,6,7,9,12 


-12000. 

13,15,16,18,19 


-12000. 

21,22,24,25,27 


-12000. 

28,30 


-12000. 

5,8,11,14,17,20 


-48000. 

23,26,29 


-48000. 

Load case 2 

Node 

x- Load (kips) 

y-Load(kips) 

1,3 


-4500. 

2, 


-18000. 

4,6,7,9,12 


-9000. 

13,15,16,18,19 


-9000. 

21,22,24,25,27 


-9000. 

28,30 


-9000. 

5,8,11,14,17,20 


-36000. 

23,26,29 


-36000. 

1 

3180. 


4 

6230. 


7 

6080. 


10 

5920. 


IS 

5740. 


16 

5540. 


19 

5300. 


22 

5000. 


25 

4610. 


28 

4010. 

Load case 3 

Node 

x-Load(kips) 

y-Load(kips) 

1,3 


-4500. 

2, 


-18000. 

4,6,7,9,12 


-9000. 

13,15,16,18,19 


-9000. 

21,22,24,25,27 


-9000. 

28,30 


-9000. 

5,8,11,14,17,20 


- 36000 . 

23,26,29 


-36000. 

1 

-3180. 


4 

-6230. 


7 

-6080. 


10 

-5920. 


13 

-5740. 


16 

-5540. 


19 

-5300. 


22 

-5000. 


25 

-4610. 


28 

-4010. 



Moment(inch • lbs) 


Moment (inch • lbs) 


Moment(inch * lbs) 


Loading Information For 40 Member Frame 
Figure 4 



STRUCTURAL WEIGHT (LBS) 


40 ELEMENT STRUCTURE 




STRUCTURAL WEIGHT (LBS) 


Load case 1 


Node 

x-Load(kips) 

y-Load(kips) 

15 


-26. 

16 


-30. 

88 


-18. 

89 


-20. 



Load case 2 

Node 

x-Load(kips) 

y-Load(kips) 

1 

2.0 


6 thru 65 by 5 

4.0 


68,75,82 

4.0 


90 thru 170 by 5 

4.0 

Load case 3 

Node 

x-Load(kips) 

y-Load(kips) 

5 

-2.0 


16,19 

-4.0 


22 thru 67 by 5 

-4.0 


74,89 

-4.0 


94 thru 174 by 5 

-4.0 



Load case 4 = Load case 1 + Load case 2 


Load case 5 — Load case 1 + Load case 3 


Momentfinch * kip) 


Moment(inch * kip) 


Moment(inch * kip) 


Loading Information For 313 Member Frame 
Figure 7 


313 ELEMENT STRUCTURE 



Weight Iteration History 
Figure 8 



